/*******************************************************************************
  ARTICLE	GAY, GOBBI, GONI (2025) "REVOLUTIONARY TRANSITIONS. INHERITANCE    
            CHANGE AND FERTILITY DECLINE" JOURNAL OF POLITICAL ECONOMY         
                                                                               
  AUTHORS	VICTOR GAY, PAULA GOBBI, MARC GONI                                 
  CONTACT	victor.gay@tse-fr.eu; paula.eugenia.gobbi@ulb.be; marc.goni@uib.no 
  VERSION	1.0 (MAY 2025)                                                     
  SOFTWARE	STATA SE 18                                                        
  LICENCE	MIT                                                                
--------------------------------------------------------------------------------

HENRY RESULTS PAPER DO FILE

This file contains the codes to reproduce the figures and tables in the paper using the Henry data.

Instructions: 
-------------
	gain access to Henry nominatif database and place enfants.sas7bdat and parents.sas7bdat in folder /1_raw_data/1_1_henri/ (see README for more details); 
	
	run 00_henry_data_prep.do once 
	
	open do-files from directory where they are placed; order matters; run whole code.
	
Contents: 
---------
	Program setup
	Figure 1: Fertility decline in France, 1700–1850.
	Figure 4: Completed fertility by cohort across inheritance systems.
	Table 1: Balancedness of pre-reform characteristics.
	Table 2: Difference-in-differences estimates, Henry database.
	Table 3: Flexible-trend difference-in-differences estimates, Henry database.
	Table 5: Heterogeneous effects by soil conditions for small versus large farms.

Date last update: May 2025; Ran using STATA 18.5
*/
********************************************************************************


********************
* 0. PROGRAM SETUP *
********************

version 18
clear all
set more off

************************
* PACKAGE DEPENDENCIES *
************************

net install tabmiss, from(https://stats.oarc.ucla.edu/stat/stata/ado/analysis/) replace
ssc install reghdfe, replace
ssc install ftools, replace
ssc install coefplot, replace
ssc install outreg2, replace

***************
* DIRECTORIES *
***************

global DAT 	= "../../3_outputs/3_1_datasets"
global TAB 	= "../../3_outputs/3_2_main/3_2_1_main_tables"
global FIG 	= "../../3_outputs/3_2_main/3_2_2_main_figures"

timer on 1

* ==============================================================================
* Figure 1: Motivation
* ------------------------------------------------------------------------------

* Henry data
use "$DAT/final-henry.dta", clear
keep if nfert>0			// sample: mothers (net)
collapse nfert, by(byear)
tsset byear
tssmooth ma mafert = nfert, window(2 1 2)
keep if byear>1650
rename byear year
replace year = year+40 // year of reference when women were aged 40
tempfile henry
save `henry', replace

* Weir Ig index and Chesnais births per 1,000
import excel "../../1_raw_data/1_18_fertility_trends/ferttrends.xlsx", sheet("Sheet1") firstrow clear
keep year NR Ig
merge 1:1 year using `henry', nogen

* Figure
gen NR0 = NR/10
sort year
drop if year==. 
keep if year>1700 
keep if year<1850
twoway 	(scatter NR0 year, m(d) mc(blue) msize(small)) ///
		(line mafert year, lc(black)) ///
	    (lfit nfert year if year< 1794, lc(black) lw(medthick) lpattern(dash)) ///
	    (lfit nfert year if year>=1794, lc(black) lw(medthick) lpattern(dash)) ///
		(scatter Ig year, yaxis(2) m(o) mc(red) msize(small)) ///
	    (pcarrowi 4.25 1793 4.25 1850 (1) " Inheritance reforms (1793)", lcolor(none) mc(none) mlabcolor(gray)) ///
		, xtitle(" ") xlabel(1700(25)1850, labs(medsmall)) xline(1793, lc(gray) lw(medthick) lp(solid)) ///
		ytitle("Births", size(medsmall)) ytitle("Ig index", axis(2) size(medsmall)) ylabel(2.5 "2.5" 3.0 "3.0" 3.5 "3.5" 4.0 "4.0" 4.5 "4.5", axis(1) labs(medsmall)) ylabel(0.5 "0.5" 0.6 "0.6" 0.7 "0.7" 0.8 "0.8" 0.9 "0.9" 1.0 "1.0", axis(2) labs(medsmall)) ///
		legend(on order(1 "Births per 100 population (INED, Chesnais 1992)" 6 "Ig index of marital fertility (Weir 1994)" 2 "Completed fertility of mothers aged 40 (Enquête Henry)") rows(3) pos(7) ring(0) size(small) symxsize(3) ) scale(1.31)
graph export "$FIG/figure1.pdf", as(pdf) name("Graph") replace


* ==============================================================================
* Figure 4: Completed fertility by cohort across inheritance systems.
* ------------------------------------------------------------------------------

* PANEL A: Raw trends
* -------------------
use "$DAT/final-henry.dta", clear
keep if byear!=. & byear>=1704 // sample
gen T2 = T if impart==1

gen a1fert = nfert if impart==1 & womexc==1
gen a0fert = nfert if impart==0 & womexc==0
gen i1fert = nfert if impart==1
gen i0fert = nfert if impart==0
gen w1fert = nfert if womexc==1
gen w0fert = nfert if womexc==0

collapse a0fert a1fert i0fert i1fert w0fert w1fert T2, by(byear)
tsset byear
label var byear "Woman's birth year"
twoway (scatter a0fert byear, m(o) mc(blue) mfc(white)) ///
	   (scatter a1fert byear, m(o) mc(red)) ///
	   (lpolyci a0fert byear if byear>=1704 & byear<=1753, lc(blue) fc(blue) fi(inten30) alc(blue%20)) ///
	   (lpolyci a1fert byear if byear>=1704 & byear<=1753, lc(red) fc(red) fi(inten30) alc(red%20)) ///
	   (lpolyci a0fert byear if byear>=1754 & byear<1804, lc(blue) fc(blue) fi(inten30) alc(blue%20)) ///
	   (lpolyci a1fert byear if byear>=1754 & byear<1804, lc(red) fc(red) fi(inten30) alc(red%20)) ///
	   (line T2 byear, yaxis(2) lc(gray%25) lw(vvthick)) ///
	   (pcarrowi 4 1754.5 4 1803 (1) " Cohorts fertile after reform (F=1)", lcolor(black) mc(black) mlabcolor(black)), ///
	   xlabel(1710(10)1800, labs(medsmall)) ///
	   ytitle("Completed fertility", size(medium)) ylabel(0(1)4, labs(medsmall) angle(horizontal)) ///
	   ytitle("Years fertile post reform", axis(2) size(medsmall) color(gray)) ///
	   ylabel(0 "0" 5 "5" 10 "10" 15 "15" 20 "20" 25 "25", labs(small) angle(horizontal) labc(gray) tlc(gray) axis(2)) yscale( lc(gray) axis(2)) ///
	   ymlabel(26 " ", notick axis(2)) ///
	   legend(on order(4 "Partible, women included" 6 "Impartible, women excluded" 12 "Years fertile post-reform (right)") cols(1) size(medsmall) pos(7) ring(0) region(lc(black))) ///
	   xline(1753.5, lcolor(black) lpattern(dash)) xtitle("Women's birth cohort", size(medium)) ///
	   aspectratio(0.5) graphregion(fcolor(white))
graph export "$FIG/figure4a.pdf", as(pdf) name("Graph") replace


* Panel B: Event study figure (R1.3c)
* -----------------------------------
use "$DAT/final-henry.dta", clear

* Sample: 
keep if byear>=1719
qui xi: reghdfe nfert T01xaffected , absorb(byear villagenum) cluster(villagenum)
keep if e(sample)==1	// regression sample

* 5-year cohort dummies
gen byearv = byear-4
gen coh = 5 * floor(byearv/5)
tab coh, gen(coh_)
foreach x of varlist coh_*{
replace  `x'= 0 if affected==0
}
replace coh_16 = coh_16+coh_17
drop coh_17 

* treatment
sort coh byear
by coh: egen Tstep = max(T)
reg Tstep coh_1-coh_16, nocons
estimates store treat0

* baseline estimate
xi: reghdfe nfert T01xaffected i.SI_MA_FE, absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE) cluster(villagenum)
estimates store base0

* event study
xi: reghdfe nfert affected coh_1-coh_6 o.coh_7 coh_8-coh_16 i.SI_MA_FE, absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE) cluster(villagenum)
estimates store affected0
coefplot (affected0, mcolor(black) ciopts(lc(black) recast(rcap)) keep(coh_*) ) ///
	(base0, mcolor(stc1) ciopts(lc(stc1) recast(rcap)) keep(T01xaffected)) ///
	(treat0, yaxis(2) noci recast(line) lc(gray%25) lwidth(vvthick) connect(stairstep)) ///
	, levels(90) name(all_grph, replace) omitted vertical yline(0, lcolor(red) lpattern(dash)) ///
	xtitle("Age at inheritance reforms (1793)") ytitle("Completed fertility") ytitle("Years fertile post-reform", axis(2) size(medium) color(gray)) ///
	xlabel(1 "75-79" 2 "70-74" 3 "65-69" 4 "60-64" 5 "50-54" 6 "45-49" 7 "40-44" 8 "35-39" 9 "30-34" 10 "25-29" 11 "20-24" 12 "15-19" 13 "10-14" 14 "5-9" 15 "0-4" 16 "<0", angle(90)) ///
	xmlabel(17 " ", notick) ///
	ylabel(-1(0.5)1, axis(1)) ///
	ylabel(#5, labc(gray) tlc(gray) axis(2)) nooffsets yscale( lc(gray) axis(2)) ///
	legend(on order(2 "Estimate by cohort" 4 "Aggregate effect" 5 "Years fertile (right)") rows(1) pos(7) ring(1) siz(medium) ) aspect(0.4)
graph export "$FIG/figure4b.pdf", as(pdf) name("Graph") replace
* ==============================================================================


* ==============================================================================
* Table 1: Balancedness of pre-reform inheritance areas.
* ------------------------------------------------------------------------------
* Data
use "$DAT/final-henry.dta", clear

* Sample
keep if byear>=1700 // born after 1700
xi: reghdfe nfert affected T01xaffected, absorb(byear) cluster(villagenum)
keep if e(sample)==1 // regression sample
keep if T==0 // keep cohorts before reform

* Variables for balance table
 // survival variables
gen wmotheralive = (SU_ME_FE==1)
gen hmotheralive = (SU_ME_MA==1)
gen wfatheralive = (SU_PE_FE==1)
gen hfatheralive = (SU_PE_MA==1)
 // accuracy variables
gen knownbdate = (TYP_FICH==11 | TYP_FICH==21 | TYP_FICH==15 | TYP_FICH==25)
gen knownudate = (TYP_FICH==11 | TYP_FICH==12 | TYP_FICH==13 | TYP_FICH==14)
* Individual-level summary stats
ssc install balancetable, replace
label var nfert "\hspace{0.7cm} Completed fertility"
label var dage "\hspace{0.7cm} Wife's age at death"
label var sdage "\hspace{0.7cm} Husband's age at death"
label var d40 "\hspace{0.7cm} Wife died before age 40"
label var wmotheralive "\hspace{0.7cm} Wife's mother alive at marriage"
label var hmotheralive "\hspace{0.7cm} Husband's mother alive at marriage"
label var wfatheralive "\hspace{0.7cm} Wife's father alive at marriage"	
label var hfatheralive "\hspace{0.7cm} Husband's father alive at marriage"
label var knownbdate "\hspace{0.7cm} Known birth year or in GPC"
label var knownudate "\hspace{0.7cm} Known union end date"
label var SI_MA_FE "\hspace{0.7cm} Literacy"
label var SI_MA_MA "\hspace{0.7cm} Literacy of husband"
global indbala dage sdage d40 wmotheralive hmotheralive wfatheralive hfatheralive knownbdate knownudate SI_MA_FE SI_MA_MA 
format nfert %4.2f
format $indbala %4.2f
balancetable (mean if affected==0) (mean if affected==1) (diff affected) ///
			 nfert $indbala using "$TAB/table1.tex", vce(cluster villagenum) ///
			 ctitles("Not reformed" "sd" "Reformed" "sd" "Difference" "se") ///
			staraux wide displayformat format(%4.2f) varlabels replace
			* missing age at death
			foreach x in "" s{
			tabmiss `x'dage if affected==0
			tabmiss `x'dage if affected==1
			}
* Municipality-level summary stats
global munbala pwheat_age15 ldensity_1793 vidxcivil ldist_eveche ldist_bailliage ldist_subdeleg ldist_recette ldist_socpol ldist_rebellion_state_1779_1789 ldist_cassini ldist_post avcalmean ruggedness sh_sandy 
collapse $munbala , by(villagenum impart womexc affected)
// 3 municipalities in alps and pyrenees: [to see type: hist ruggedness, bin(39) freq ylabel(0(1)8) ytitle("Terrain ruggedness (in 100s of meters)")]
replace ruggedness = . if ruggedness>3
label var pwheat_age15 "\hspace{0.7cm} Wheat price (log)"
label var ldensity_1793 "\hspace{0.7cm} Population density"
label var vidxcivil "\hspace{0.7cm} Religiosity index"
label var ldist_eveche "\hspace{0.7cm} Distance to religious centre (Eveche)"
label var ldist_bailliage "\hspace{0.7cm} Distance to judicial district seat"
label var ldist_subdeleg "\hspace{0.7cm} Distance to territorial admin. (subdelegation)"
label var ldist_recette "\hspace{0.7cm} Distance to tax centre (Recette)"
label var ldist_socpol "\hspace{0.7cm} Distance to political society"
label var ldist_rebellion_state_1779_1789 "\hspace{0.7cm} Distance to rebellion in 1779-89"
label var ldist_cassini "\hspace{0.7cm} Distance to Cassini road in 1750-90"
label var ldist_post "\hspace{0.7cm} Distance to horse post in 1790"
label var avcalmean "\hspace{0.7cm} Average caloric suitability of land"
label var ruggedness "\hspace{0.7cm} Terrain ruggedness (in 100s of meters)"
label var sh_sandy "\hspace{0.7cm} Share sandy soils"

format $munbala %4.2f
balancetable (mean if affected==0) (mean if affected==1) (diff affected) ///
			 $munbala using "$TAB/table1.tex", vce(robust) ///
			 ctitles("Not reformed" "sd" "Reformed" "sd" "Difference" "se") ///
			 staraux wide displayformat format(%4.2f) varlabels append
* ==============================================================================


* ==============================================================================
* Table 2: Difference-in-differences estimates, Henry database.
* ------------------------------------------------------------------------------
use "$DAT/final-henry.dta", clear

* Sample
keep if byear>=1700 // born after 1700
qui xi: reghdfe nfert T01xaffected, absorb(byear) cluster(villagenum)
keep if e(sample)==1 // reg non-missing values

* Regressions
label var T01xaffected "Reform $\times$ Fertile"
local x1 "Dep. var.: completed fertility"
local x2 "Dep. var.: completed fertility of mothers"
local x3 "Dep. var. = 1 if childless"
local x4 "Dep. var. = age at marriage"
local i=1
local panels="A B C"
foreach x in nfert intmarg chdls{
global depvar `x'
local pan `: word `i' of `panels''

// (1) no controls
xi: reghdfe $depvar T01xaffected, absorb(byear villagenum) cluster(villagenum)
outreg2 using "$TAB/table2-panel`pan'.tex", replace ///
	keep(T01xaffected) nocons ///
	title("`x`i''") ctitle("`x'") ///
	addtext(Cohort FE, Y, Municipality FE, Y, Literacy, ., Accuracy of Henry form FE, ., Father alive at marriage, ., Mother alive at marriage, .) ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long
// (2) + literacy
xi: reghdfe $depvar T01xaffected i.SI_MA_FE , absorb(byear villagenum) cluster(villagenum)
outreg2 using "$TAB/table2-panel`pan'.tex",  ///
	keep(T01xaffected) nocons ///
	title("`x`i''") ctitle("`x'") ///
	addtext(Cohort FE, Y, Municipality FE, Y, Literacy, Y, Accuracy of Henry form FE, ., Father alive at marriage, ., Mother alive at marriage, .) ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long
// (3) + FE for dates' accuracy (type fiche)
xi: reghdfe $depvar T01xaffected i.SI_MA_FE, absorb(byear villagenum TYP_FICH) cluster(villagenum)
outreg2 using "$TAB/table2-panel`pan'.tex",  ///
	keep(T01xaffected) nocons ///
	title("`x`i''") ctitle("`x'") ///
	addtext(Cohort FE, Y, Municipality FE, Y, Literacy, Y, Accuracy of Henry form FE, Y, Father alive at marriage, ., Mother alive at marriage, .) ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long
// (4) + parental survival at marriage
xi: reghdfe $depvar T01xaffected i.SI_MA_FE, absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE) cluster(villagenum)
outreg2 using "$TAB/table2-panel`pan'.tex",  ///
	keep(T01xaffected) nocons ///
	title("`x`i''") ctitle("`x'") ///
	addtext(Cohort FE, Y, Municipality FE, Y, Literacy, Y, Accuracy of Henry form FE, Y, Father alive at marriage, Y, Mother alive at marriage, Y) ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long

local i = `i'+1
}
* ==============================================================================


* ==============================================================================
* Table 3: Flexible-trend difference-in-differences estimates, Henry database.
* ------------------------------------------------------------------------------
use "$DAT/final-henry.dta", clear

* Sample: born after 1700
keep if byear>=1700

* Regressions
label var T01xaffected "Reform $\times$ Fertile"
local x1 "Dep. var.: completed fertility"
local x2 "Dep. var.: completed fertility of mothers"
local x3 "Dep. var. = 1 if childless"
global flextrends prewheat ldensity_1793 vidxcivil ldist_eveche ldist_socpol ldist_rebellion_state_1779_1789 ldist_bailliage ldist_subdeleg ldist_recette ldist_cassini ldist_post
local i=1

* Completed fertility
// + wheat prices
reghdfe nfert T01xaffected i.SI_MA_FE logpwheat, absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.(prewheat ldensity_1793)#byear) cluster(villagenum)
keep if e(sample)==1
outreg2 using "$TAB/table3.tex", ///
	keep(T01xaffected) nocons ///
	title("flexible trends") ctitle("CF") ///
	addtext(Cohort and municipality FE, Y, Individual-level controls, Y, Local wheat price in decade logs, Y, Initial wheat price $\times$ Cohort FE, Y, Initial log pop density $\times$ Cohort FE, Y, Religiosity index $\times$ Cohort FE, ., Distance religious center $\times$ Cohort FE, ., Distance political society $\times$ Cohort FE, ., Distance rebellion 1779-89 $\times$ Cohort FE, ., Distance legal center $\times$ Cohort FE, ., Distance fiscal center $\times$ Cohort FE, ., Distance administrative center $\times$ Cohort FE, ., Distance Cassini road $\times$ Cohort FE, ., Distance horse post $\times$ Cohort FE, .) ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long replace
	
// + flexible-trends by religiosity
reghdfe nfert T01xaffected i.SI_MA_FE logpwheat, absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.(prewheat ldensity_1793 vidxcivil ldist_eveche)#byear) cluster(villagenum)
outreg2 using "$TAB/table3.tex", ///
	keep(T01xaffected) nocons ///
	title("flexible trends") ctitle("CF") ///
	addtext(Cohort and municipality FE, Y, Individual-level controls, Y, Local wheat price in decade logs, Y, Initial wheat price $\times$ Cohort FE, Y, Initial log pop density $\times$ Cohort FE, Y, Religiosity index $\times$ Cohort FE, Y, Distance religious center $\times$ Cohort FE, Y, Distance political society $\times$ Cohort FE, ., Distance rebellion 1779-89 $\times$ Cohort FE, ., Distance legal center $\times$ Cohort FE, ., Distance fiscal center $\times$ Cohort FE, ., Distance administrative center $\times$ Cohort FE, ., Distance Cassini road $\times$ Cohort FE, ., Distance horse post $\times$ Cohort FE, .) ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long
	
// + flexible-trends by political factors + economic-geography
reghdfe nfert T01xaffected i.SI_MA_FE logpwheat, absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.(prewheat ldensity_1793 vidxcivil ldist_eveche ldist_socpol ldist_rebellion_state_1779_1789 ldist_bailliage ldist_subdeleg ldist_recette ldist_cassini ldist_post)#byear) cluster(villagenum)
outreg2 using "$TAB/table3.tex", ///
	keep(T01xaffected) nocons ///
	title("flexible trends") ctitle("CF") ///
	addtext(Cohort and municipality FE, Y, Individual-level controls, Y, Local wheat price in decade logs, Y, Initial wheat price $\times$ Cohort FE, Y, Initial log pop density $\times$ Cohort FE, Y, Religiosity index $\times$ Cohort FE, Y, Distance religious center $\times$ Cohort FE, Y, Distance political society $\times$ Cohort FE, Y, Distance rebellion 1779-89 $\times$ Cohort FE, Y, Distance legal center $\times$ Cohort FE, Y, Distance fiscal center $\times$ Cohort FE, Y, Distance administrative center $\times$ Cohort FE, Y, Distance Cassini road $\times$ Cohort FE, Y, Distance horse post $\times$ Cohort FE, Y) ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long

* Completed fertility of mothers
// + all flexible-trends
reghdfe intmarg T01xaffected i.SI_MA_FE logpwheat, absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.($flextrends )#byear) cluster(villagenum)
outreg2 using "$TAB/table3.tex", ///
	keep(T01xaffected) nocons ///
	title("flexible trends") ctitle("CF mothers") ///
	addtext(Cohort and municipality FE, Y, Individual-level controls, Y, Local wheat price in decade logs, Y, Initial wheat price $\times$ Cohort FE, Y, Initial log pop density $\times$ Cohort FE, Y, Religiosity index $\times$ Cohort FE, Y, Distance religious center $\times$ Cohort FE, Y, Distance political society $\times$ Cohort FE, Y, Distance rebellion 1779-89 $\times$ Cohort FE, Y, Distance legal center $\times$ Cohort FE, Y, Distance fiscal center $\times$ Cohort FE, Y, Distance administrative center $\times$ Cohort FE, Y, Distance Cassini road $\times$ Cohort FE, Y, Distance horse post $\times$ Cohort FE, Y) ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long

* Childlessness
// + all flexible-trends
reghdfe chdls T01xaffected i.SI_MA_FE logpwheat, absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.($flextrends )#byear) cluster(villagenum)
outreg2 using "$TAB/table3.tex", ///
	keep(T01xaffected) nocons ///
	title("flexible trends") ctitle("Childless") ///
	addtext(Cohort and municipality FE, Y, Individual-level controls, Y, Local wheat price in decade logs, Y, Initial wheat price $\times$ Cohort FE, Y, Initial log pop density $\times$ Cohort FE, Y, Religiosity index $\times$ Cohort FE, Y, Distance religious center $\times$ Cohort FE, Y, Distance political society $\times$ Cohort FE, Y, Distance rebellion 1779-89 $\times$ Cohort FE, Y, Distance legal center $\times$ Cohort FE, Y, Distance fiscal center $\times$ Cohort FE, Y, Distance administrative center $\times$ Cohort FE, Y, Distance Cassini road $\times$ Cohort FE, Y, Distance horse post $\times$ Cohort FE, Y) ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long
	
* Age at marriage
// + all flexible-trends
reghdfe mage T01xaffected i.SI_MA_FE logpwheat, absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.($flextrends )#byear) cluster(villagenum)
outreg2 using "$TAB/table3.tex", ///
	keep(T01xaffected) nocons ///
	title("flexible trends") ctitle("Age at marriage") ///
	addtext(Cohort and municipality FE, Y, Individual-level controls, Y, Local wheat price in decade logs, Y, Initial wheat price $\times$ Cohort FE, Y, Initial log pop density $\times$ Cohort FE, Y, Religiosity index $\times$ Cohort FE, Y, Distance religious center $\times$ Cohort FE, Y, Distance political society $\times$ Cohort FE, Y, Distance rebellion 1779-89 $\times$ Cohort FE, Y, Distance legal center $\times$ Cohort FE, Y, Distance fiscal center $\times$ Cohort FE, Y, Distance administrative center $\times$ Cohort FE, Y, Distance Cassini road $\times$ Cohort FE, Y, Distance horse post $\times$ Cohort FE, Y) ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long
* ==============================================================================


* ==============================================================================
* Table 5. HETEROGENEITY BY SOIL TEXTURE
* ------------------------------------------------------------------------------
use "$DAT/final-henry.dta", clear

set seed 03122024

global flextrends prewheat ldensity_1793 vidxcivil ldist_eveche ldist_socpol ldist_rebellion_state_1779_1789 ldist_bailliage ldist_subdeleg ldist_recette ldist_cassini ldist_post

* sample
keep if byear>=1700 // born after 1700
qui reghdfe nfert T01xaffected i.SI_MA_FE logpwheat, nocons absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.($flextrends )#byear) cluster(villagenum)
keep if e(sample)==1 // regression sample

* Heterogeneous effects DD regressions
gen T01xaffectedXsclc1 = .
gen T01xaffectedXsclc0 = .
label var T01xaffectedXsclc1 "Reform X Fertile (soil conditions for large farms)"
label var T01xaffectedXsclc0 "Reform X Fertile (soil conditions for small farms)"
global treatvars T01xaffectedXsclc1 T01xaffectedXsclc0

// Based on soil texture
replace T01xaffectedXsclc1 = T01xaffected*(text1_cat==1)
replace T01xaffectedXsclc0 = T01xaffected*(text1_cat!=1)
* col (1)
qui reghdfe nfert $treatvars i.SI_MA_FE logpwheat, nocons absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.($flextrends )#byear) cluster(villagenum)
local nobs = `e(N)'
local diff = round(_b[T01xaffectedXsclc1] - _b[T01xaffectedXsclc0],0.001)
test _b[T01xaffectedXsclc1] - _b[T01xaffectedXsclc0] = 0
local pval = round(`r(p)',0.001)
outreg2 using "$TAB/table5.tex", ///
	keep($treatvars ) nocons noobs nor2 ///
	addtext(Differential effect B1 B2, `diff', p value difference, `pval', Cohort FE, Y, Municipality FE, Y, Individual-level controls, Y, Flexible trends, Y, Soil texture X cohort FE, N, Ruggedness X cohort FE, N, Soil conditions based on, Soil texture, N, `nobs') ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long replace
* col (2)	
qui reghdfe nfert $treatvars i.SI_MA_FE logpwheat, nocons absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.($flextrends text1 ruggedness)#byear) cluster(villagenum)
local nobs = `e(N)'
local diff = round(_b[T01xaffectedXsclc1] - _b[T01xaffectedXsclc0],0.001)
test _b[T01xaffectedXsclc1] - _b[T01xaffectedXsclc0] = 0
local pval = round(`r(p)',0.001)
outreg2 using "$TAB/table5.tex", ///
	keep($treatvars ) nocons noobs nor2 ///
	addtext(Differential effect B1 B2, `diff', p value difference, `pval', Cohort FE, Y, Municipality FE, Y, Individual-level controls, Y, Flexible trends, Y, Soil texture X cohort FE, Y, Ruggedness X cohort FE, Y, Soil conditions based on, Soil texture, N, `nobs') ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long

// Based on terrain ruggedness
qui bys insee_com: gen ruggedness_munlevel = ruggedness if _n==1
qui summ ruggedness_munlevel, d
qui gen iruggedness = (ruggedness<=`r(p50)')
replace T01xaffectedXsclc1 = T01xaffected*(iruggedness)
replace T01xaffectedXsclc0 = T01xaffected*(1-iruggedness)
* col (3)
qui reghdfe nfert $treatvars i.SI_MA_FE logpwheat, nocons absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.($flextrends )#byear) cluster(villagenum)
local nobs = `e(N)'
local diff = round(_b[T01xaffectedXsclc1] - _b[T01xaffectedXsclc0],0.001)
test _b[T01xaffectedXsclc1] - _b[T01xaffectedXsclc0] = 0
local pval = round(`r(p)',0.001)
outreg2 using "$TAB/table5.tex", ///
	keep($treatvars ) nocons noobs nor2 ///
	addtext(Differential effect B1 B2, `diff', p value difference, `pval', Cohort FE, Y, Municipality FE, Y, Individual-level controls, Y, Flexible trends, Y, Soil texture X cohort FE, N, Ruggedness X cohort FE, N, Soil conditions based on, Ruggedness, N, `nobs') ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long
* col (4)
qui reghdfe nfert $treatvars i.SI_MA_FE logpwheat, nocons absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.($flextrends text1 ruggedness)#byear) cluster(villagenum)
local nobs = `e(N)'
local diff = round(_b[T01xaffectedXsclc1] - _b[T01xaffectedXsclc0],0.001)
test _b[T01xaffectedXsclc1] - _b[T01xaffectedXsclc0] = 0
local pval = round(`r(p)',0.001)
outreg2 using "$TAB/table5.tex", ///
	keep($treatvars ) nocons noobs nor2 ///
	addtext(Differential effect B1 B2, `diff', p value difference, `pval', Cohort FE, Y, Municipality FE, Y, Individual-level controls, Y, Flexible trends, Y, Soil texture X cohort FE, Y, Ruggedness X cohort FE, Y, Soil conditions based on, Ruggedness, N, `nobs') ///
	addstat(N clusters, `e(N_clust)') adjr2 nonotes addnote(SE clustered by village, *p<.05; **p<.01; ***p<.001) ///
	label dec(3) adec(0) tex(fragment) long
	
* Statement in ft 50 "The sample includes three villages in the Alps and Pyrenees with extreme ruggedness. When dropping these, we obtain coefficients of −0.56 to −0.74 for rugged terrains and of −0.00 to 0.04 for flat terrains, along with a p-value on the difference of 0.004 to 0.07."	
reghdfe nfert $treatvars i.SI_MA_FE logpwheat if ruggedness<3, nocons absorb(byear villagenum TYP_FICH SU_PE_FE SU_ME_FE c.($flextrends text1 ruggedness)#byear) cluster(villagenum)
test _b[T01xaffectedXsclc1] - _b[T01xaffectedXsclc0] = 0
* ==============================================================================

timer off 1 /* 14 seconds */
timer list	